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ABSTRACT 

A method for inverting the statistical star counts equation, including proper motions, is pre- 
sented; in order to break the degeneracy in that equation it uses the supplementary constraints 
required by dynamical consistency. The inversion gives access to both the kinematics and the 
luminosity function of each population in three regimes: the singular ellipsoid, the constant ra- 
tio Schwarzschild ellipsoid plane parallel models and the epicyclic model. This more realistic 
model is taylored to account for local neighbourhood density and velocity distribution. 

The first model is fully investigated both analytically and via means of a non-parametric 
inversion technique, while the second model is shown to be formally its equivalent. The effect 
of noise and incompleteness in apparent magnitude is investigated. The third model is inves- 
tigated via a 5D+2D non-parametric inversion technique where positivity of the underlying 
luminosity function is explicitely accounted for. 

It is argued that its future application to data such as the Tycho catalogue (and in the 
upcoming satellite GAIA) could lead - provided the vertical potential, and/or the asymmet- 
ric drift or wq are known - to a non-parametric determination of the local neighbourhood 
luminosity function without any reference to stellar evolution tracks. It should also yield the 
proportion of stars for each kinematic component and a kinematic diagnostic to split the thin 
disk from the thick disk or the halo. 

Key words: methods: data analysis stars: Hertzsprung-Russell, luminosity function Galaxy: 
kinematics and dynamics, structure, stellar content 



1 INTRODUCTION 

Most of our knowledge of the global structure of the Galaxy relies 
on the comparison of magnitude and colour star counts in different 
Galactic directions. Star counts alone do not allow us to solve the 
dilemma that a star of a given apparent magnitude can be either 
intrinsically faint and close by, or bright and distant. This problem 
may be addressed statistically by using the century-old equation of 
stellar statistics (von Seeliger 1898): 

p oo 

A x (mJ,b) = / $x(M)p(r,£,b)r 2 dr, (1) 
Jo 

where A\(m, £, b) dm did( sin b) is the number of stars 
which have an apparent magnitude in the range [m, m + dm] , 
&\(M) is the luminosity function, which depends on the intrin- 
sic magnitude, M, and the colour band A, while p(r,£,b) is the 
density at radius r (within dr) along the line of sight in the direc- 
tion given by the Galactic longitudes and latitudes (£, b) (within the 
solid angle d£ cos(6)db). 

This equation cannot be solved or inverted (i.e. by determining 
both the stellar LF and the density law) except for a few simplified 
cases. For instance, with a "homogeneous" stellar sample for which 



the absolute magnitudes of stars or more precisely their luminosity 
functions are known, the density law along the line of sight can 
be recovered. A classical numerical technique (Mihalas & Binney 
1981) has been proposed - the Bok diagram (1937) - while more 
rigorous treatments are required for small samples to stabilize the 
inversion so as to produce smooth solutions (Binney & Merrifield 
1998). The converse situation is the determination of the luminosity 
function assuming a known density law (see, for instance, recent 
studies of the faint end of the disk or halo main sequence based on 
deep star counts (Reid et al. 1996, Gould et al. 1998). 

A simple approach, developed largely in the eighties, was to 
integrate Eq. ([j]) assuming some prior information concerning the 
stellar populations (see, for instance, Pritchet 1983, Bahcall et al., 
1983, Buser 1985, Robin & Creze 1986). A frequent assumption is, 
for instance, to assume that the halo stars have the same luminos- 
ity function as some low metallicity globular clusters. Another ap- 
proach consists in building a stellar luminosity function from stellar 
evolution tracks and isochrones of various ages. This has been used 
to put constraints on the Galactic disk star formation rate (Haywood 
et al 1997ab). 
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Stronger a priori constraints may also be derived by requir- 
ing dynamical consistency, since the vertical kinematics of stars is 
related to the flattening of stellar disks or spheroidal components. 

Since star counts alone, A\(m, £, b), are not sufficient to con- 
strain uniquely Galactic stellar population models, it is expected 
that two (or more) distinct models will reproduce the same appar- 
ent star counts. However, this is not a real worry, since it is likely 
that adding some relevant extra a priori information must help to 
lift partially the degeneracy of the models. 

In this paper, it is shown that the degeneracy is lifted alto- 
gether when we consider, in addition to the star counts in apparent 
magnitude, the proper motions, \ii and fit,- For a relatively general 
dynamically consistent model (stationary, axisymmetric and fixed 
kinematic radial gradients), the statistical equation counts may be 
formally inverted, giving access to both the vertical density law of 
each stellar population and their luminosity functions. This is de- 
veloped in section ^ where we show how the vertical motions are 
related to the thickness of stellar components. The remaining de- 
generacy occurs only for a quadratic vertical potential. Otherwise 
— when the vertical component of the potential is known — the de- 
partures from quadratic behaviour define a characteristic scale that 
allows us to transform statistically the magnitudes into distances 
and proper motions into velocities. Similarly, the asymmetric drift 
and/or the vertical velocity component of the Sun provide a natural 
scale in energy, leading to the same inversion procedure. 

For ideal star counts (infinitely deep and for an infinite number 
of stars), the inversion gives exactly the proportion of stars in each 
kinematic component, providing a direct diagnostic to split the thin 
disk from the thick disk or the halo, and its luminosity function 
$a (M) is recovered for each kinematic stellar component. This is 
a direct consequence of the supplementary constraints introduced 
by the requirement for dynamical consistency. 

Section ^ presents the generalized stellar statistic equation 
which accounts for proper motions, and demonstrates the unique- 
ness of the inversion for two families of plane parallel distribution 
functions: the singular velocity ellipsoid (subsection 2. 1 ) and a con- 
stant ratio velocity ellipsoid (subsection 2.2 ) while subsection 23 
presents a basic description of the epicyclic model. Section p il- 
lustrates the inversion procedure on a fictitious superposition of 4 
kinematically decoupled populations with distinct main sequence 
turn-off magnitudes for the constant ratio velocity ellipsoid and the 
epicyclic models. The next section discusses the effects of trunca- 
tion in apparent magnitude (i.e. completude of the catalog) in the 
recovered LF as well as noise in the measurements. Finally, the 
last section discuss the applicability of the method to the Tycho-2 
catalogue and to external clusters, and concludes. 
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2 DERIVATION 

The number of stars, AN, which have an apparent luminosity in the range [L, L + AL[ in the solid angle defined by the Galactic longitudes 
and latitudes (I, b) (within d^d(sin 6)), with proper motions \n and fit (within Afib and Afie ) is given by 

AN = A x (L,fn,fi b ;l,b)Afi e Afi b A£cosbAbAL=< / / $* X [L O ,0] I fp(r, u)Au r ) r 4 drd/3 > Afi t Afi b A£ cos bAbAL , (2) 



where we have introduced the luminosity function per unit bandwidth, $J [Lo,/S], which is here taken to be a function of the absolute 
luminosity, Lo, and of a continuous kinematic index, /3. The variables r, u are the vector position and velocity coordinates (u r ,up,u b ) 
in phase space relative to the Local Standard of Rest, while R and V are those relative to the Galactic centre. The relationship between 
A\(L, fii,fib', I, b) and [Lq, /3] involves a double summation over /3, and distance, r, along the line of sight. Here fp(r, u) represents the 
component of the distribution function of the assumed stationary axisymmetric equilibrium, i.e. 

/(p,u)= / Mr, u) A/3 (3) 
J o 

where / is decomposed over the basis of isothermal solutions fp of the Boltzmann equation for the assumed known potential %p. Eq. 
corresponds to a decomposition over isothermal populations of different kinematic temperatures, a 2 = 1/(3. Apart from restriction, the shape 
of the distribution f{E z ) could be anything. Note that Eq. (^) is a direct generalization of Eq. ([[]) since 

p 8 (r,£,b,H£,fj, b ) = r 2 Jff3(r,u)Au r 

is by definition the density of stars (belonging to population j3) which are at position r = (r, £, b) within drd^d(sin b), with proper 
motion fib (within Afi b ) and /w(within Afie). The extra summation on f3 which arises in Eq. (Q) accounts for the fact that stars in the local 
neighbourhood come from a superposition of different kinematic populations which, as is shown later, can be disentangled. Note that <&* x is 
defined here per unit absolute luminosity, Lo, and therefore 

$ A (M(L )) = 21 ° S 5 (10) Lo J $t[L ,l3]AP where M(L ) = -5_i_iog (^^j + M Q . 

Since there is no convolution on A (which is mute) it will be omitted from now on in the derivation. In section Q B-V colours are 
reintroduced to demonstrate the inversion for a fictitious HR diagram. We shall also drop the * superscript but will keep in mind that the 
luminosity function is expressed as a function of the absolute luminosity, Lq. 

This paper is concerned with the inversion of Eq. We proceed in three steps; first a simplistic Ansatz for the distribution function is 
assumed (corresponding to a stratification in height of uniform disks with a pin-like singular velocity ellipsoid) leading to a proof that, in this 
context, Eq. (^) has a well-defined unique solution which can be made formally explicit. A more realistic model is then presented accounting 
for the measured anisotropy of the velocity ellipsoid. It is shown that, in the direction of the Galactic centre, and if the velocity dispersions 
ratios are constant for all populations, this model is formally invertible following the same route. Away from the Galactic centre direction, the 
Sun's velocities are also accounted for to recover statistically distances via another inversion procedure related to secular parallaxes. Finally, 
we illustrate the inversion on a fully 7 dimensional epicyclic model. The detailed investigation of this model is postponed to a companion 
paper, (Siebert, Pichon, & Bienayme in preparation). 



2.1 A toy model: parallel sheet model with singular velocity ellipsoid 

Let us assume here a sheet-like model for the distribution function of kinematic temperature /3: 

Mr, u) = ^ ex P S(vr)S(v^) , (4) 

which corresponds to a stratification in height with a pin-like singular velocity ellipsoid which is aligned with the rotation axis of the Galaxy. 
Calling fib = Ub /r, the energy reads in terms of the heliocentric coordinates 

= f + M*) = + «sin 2 (6)r 2 + x (r sin(b)) , (5) 

where the harmonic component of the z potential (az 2 ) was made explicit while leaving unspecified the non harmonic-residual, \- 
Putting Eq. (§) into Eq. (|) leads to 

A[b,n b ,L) = J J exp ^«r 2 sin 2 (b) -/3r 2 2co M j (b) - fo(r sin(6))^ r 3 ArAf3 , (6) 

given the relationship Lo = Lr 2 relating apparent and absolute luminosities. Introducing £ = L 1//2 r , 
sin 2 (6) u? . sin(6) 

x = a ^ + 2i^m> and y =iM- (7) 

Eq. da) then reads 
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L 2 cos(b)A[b, fj,b, L] 



: $ [C 2 ,/3] exp [~(3( 2 x - Px(Cyj\ ( 3 d(d(3. 



(8) 



2.1.1 harmonic degeneracy 

Suppose for now that the z-potential is purely harmonic, so that \ is identically null. Calling s — /3£ 2 , the inner integral over £ in Eq. 
can be rewritten as an integral over /3 and s 



^$[C 2 ,/3]exp[-/3< 2 x]< 3 dCd/3 1 



2V2tt 



$ [s/0, 13] /3" 3/2 d/3 ) exp [sx] sds 



(9) 



Eq. (JpJ) shows that for a purely harmonic potential the mixture of populations (integrated over /3) is recovered from A[b, fib, L] which is 
effectively a function of x only (given by Eq. (Q)). In this instance, the inversion does not allow us to disentangle the different kinematic 
populations. In physical terms there is a degeneracy between the distance, luminosity and proper motion. In contrast, when the data set 
extends far enough to probe the an-harmonic part of the potential, we now demonstrate that Eq. (^) has formally a unique exact solution, 
before exploring non-parametric means of inverting it in a more general framework. 



2.1.2 Uniqueness ? 

Let us assume that not too far from the Galactic plane, x( 2 ) is we ll approximated by x( 2 ) = l zl ' ■> so that Eq. (|j) becomes 

T 



L 2 cos(b)A[b,fj,b,L] 
Calling 



<J> [C 2 , P] exp (-pifx - /3 7 C V) C 3 d<d/3 . 



$i [U, B] 



1 



$ [exp(2C/), exp(B)] exp(4[/ + 3/2B - c([2 + v]U + 25)), K (z) = exp(cz - exp(z)) , 



(10) 

(11) 
(12) 

(13) 
(14) 

The positive scalar c is left to our discretion and can be chosen so as to yield a narrow kernel, Ko (in practice c should be close to one). Since 
r runs from zero to infinity and so does /3, the integration over B and Z will run from — oo to oo. Similarly X and Y span ] — oo, oo[ as b 
goes from zero to 7r/2. Let 



and Ai [X, Y] = L 2 cos bA[b, fj, b , L] exp [c{X + Y)] 
where 

sin 2 (6) 

log 



+ 



2Lcos 2 (fe) 



Y = loe 



B = log(/3), Z = log(C), ^ = log(x; 
Eq. (|ic|) becomes 

Ai^/Xi.L] = ArlX.r] = / [ 3>i[Z,B]K (B + 2Z + X)K (B + vZ + Y)dZdB 



7 sin" (6) 



W = -(B + vZ), m = -(B + 2Z), 
Eq. (0) reads 



(17) 



(15) 



Ai[X,r] = \v-2\~ 1 J / $i[a7,w]Ko(Jf-«7)iii:o(y-iu)dti7d«). (16) 
The unique solution of Eq. (Jlfj) reads formally 
$i [uj,w] = \v-2\ FT 
where 

Ax \h^j , /c^] : 
while 

FT^ 1 (/(fcx, fe a )) = y y ex P[ — — ifej/H / (k x ,ky) dkxdk y . 

Both Fourier transforms are well-defined given the span of zu, w and X, Y. Approximating both Ko and A\ by a Gaussian of width 
respectively 1/Ik and Eq. ( fl7| ) shows that $i will be a Gaussian of width 1/(£n — £k) 2 - 

This procedure is therefore a true deconvolution: the luminosity function Q\(L,P) is effectively recovered at arbitrary resolution (in 
effect fixed by the signal to noise ratio of the data). In practi ce, Eq. ( |r7| ) is impractical for noisy finite data sets, so we shall investigate 
non-parametric regularised solutions to Eq. M) in section 3.1.1 



A\ [fcro j fcw\ 

K (k m )K (k u 

exp[+i(k w X + k za Y)}Ai[X,Y]dXdY and K [k] = J exp [+i QeX)] K (X) dX , 
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There is a natural scale Bo = ( a /7) given by the break in the potential which provides us with a means to lift the degeneracy 

between faint close stars moving slowly and bright stars moving faster farther out. This scale reflects the fact that statistically the dynamics 
(i.e. the velocities) gives us a precise indication of distances in units of to. We can therefore reassign a posteriori distances to stars in the 
statistical sense and deconvolve the colour magnitude diagram. Fig. (|l|) graphically demonstrates the requirement to access the break radius 
of the potential in order to derive statistical distances to the stars. It shows sections of increasing apparent magnitude in the b, fib plane for a 
two- temperature model and for a one-temperature model (corresponding to a unique absolute luminosity). The observed star counts enable 
us to distinguish between the one and two-temperature models especially at the faint end for significantly non-zero 7. 

Turning back to Eq. M, it remains that for more general \ it can still be inverted in the least squares sense, but this involves a less 
symmetric kernel, K\ (x, y]u, (3), whose functional form depends explicitly on \'- 

K 1 (x,y\u,/3) = -y^exp [-(3u 2 x - Px{uy)] u 3 . 

The inversion procedure which will be described in section Bl still applies to such kernels. 



2.2 A Schwarzschild Model: accounting for the local velocity ellipsoid anisotropy 

Let us now move to more realistic models with a fully triaxial Schwarzschild Ellipsoid. Its distribution function is given in terms of the 
kinematic inverse dispersions /3r, (5$ and /3 Z by: 



Mr, u) = y exp (-09.g, + p R E R + , (18) 

where 

E z = i« 2 2 + 4> z (z), E R =-v 2 R and = - (vj, - v^f . (19) 

Here Vj, measures the mean azimuthal velocity in the local neighbourhood (which is assumed not to depend on /3), V = (vr, v$, v z ) are 
respectively the radial, azimuthal and vertical velocities of a given star measured in a direct cylindrical system of coordinates centred at the 
Galactic centre. These velocities are given as a function of the velocities measured in the frame of the sun by 

= sin(fc) sin(£)ub — tq cos (b) sm(£)u r — r© cos(£)ut + r cos (b) (lit — sin(£)«0) + (tq + r cos (b) cos(£)) vq) , (20) 

vr — -j- {(r cos(6) — tq cos(£)) sin(6) ut — vq s'm(£) U( — cos(&) (r cos(6) — r© cos(£)) u r + 

r© uq — r cos(b) cos{£) uq + r cos(6) sin(^) vq } , (21) 
Vz — sin(b) u r + cos(6) Ub + wq , (22) 

where 



R — yr| - 2r©r cos(fe) cos(^) + r 2 cos(6) 2 and z = r sin(6) . (23) 

R measures the projected distance (in the meridional plane) to the Galactic centre, while 2 is the height of the star. Here uq,vq, wq and r© 
are respectively the components of the Sun's velocity and its distance to the Galactic centre. The argument of the exponential in Eq. ( |l8| ) is a 
quadratic function in u r via Eqs. (H)-(^2]), so the integration over that unknown velocity component is straightforward. 

In short, we show in Appendix^ that Eq. (Q) has solutions for families of distribution obeying (JT^). Those solution are unique and can 
be made explicit for a number of particular cases which are discussed there. They are shown to be formally equivalent to those found for 
Eq. (^). For instance, at large distance from the Galactic centre (r© — ► 00) Eq. (^) along the plane fit — can be re-casted into 

L 2 cos(b)A 2 [b,£,fi b = 0,L] = ff [u,P] exp (-/3« 2 x 3 - f3z 2 ) u 3 dud/3 , with x 3 = a Sm ^ , (24) 



and 



_ (w Q cos b-[v Q -v<p] sin b sin if _ _ /3* 

22 ~ 2cos 2 (fe) + 2sin 2 (6)(^cos 2 (Q+^sin 2 (^)) ' Pr ~ ^ ' 



which is of the form described in section 



2.1.2) with v = 0, £3 replacing x and Z2 replacing y. We the exeption of the special cases also 



described in Appendix the solution can be found via \ 2 minimisation as shown below in section ^| 



2.3 Epicyclic Model : Accounting for density gradients 

The above models do not account for any density or velocity dispersion gradients, which is a serious practical shortcoming. Let us therefore 
construct an epicyclic model for which the radial variation of the potential and the kinematic properties of the Galaxy are accounted for. 

A distribution function solution of Boltzmann equation with two integrals of motion (energy and angular momentum) can be writen 
according to Shu (1969) as 





Figure 1. In each panel: Isocontour of A*(b, pt,) (defined by Eq. (|48|)) in the b, pt, plane (b ranging from —pi/2 to tt/2 and p^ from —1 to 1): sections of 
increasing apparent magnitude (from left to right and top to bottom) Top Left: two-temperature models (log j3 = — 2andlog/3 = 2,log(Lo) = 0) TopRight: 
same as top left, but for a unique temperature (/3 = 1) model. The observed star counts enable us to distinguish between the one and two-temperature models 
especially at the faint end (top left section) for significantly non zero 7 (= 1). Bottom Left: Shows that even the faint end (top left section) of the observed 
star counts are barely distinguishable from the two-temperature model (Bottom Right:) for small 7 (= 1/10) . This demonstrates graphically the requirement 
to access the break radius, £q oc I/7, of the potential to derive statistical distances to the stars. 



where O is the Heaviside function while 

Rq — R c 



fi = , p D = p exp " , with R c = H—i R-+ 1 V Q a+1 , (27) 

y/2 a + 2 \ tip J 

a being the slope of the rotation curve, Q the angular velocity, k the epicyclic frequency, po the density, R c the radius of the circular orbit 
of angular momentum H, o~% and a\ the square of the radial and vertical velocity dispersion and j3 the kinematic index 



22 I 2J?0 — 2R C \ 2 2 / 2Rq — 2R C \ „ a + 1 D -^qrr T r -^rr , 

^js = cTfi exp I - I , cr z = cr Z0 exp I - I , £ c = R Q + V Q + . (28) 



Here po,tt, k, o~r, a z and _E C are known functions of momentum H given by 
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Figure 2. Left; Aitof projection of the normalised density distribution for the epicyclic Shu model. Right: Distribution of the maximum of the proper motion 
along the galactic longitude (in "/yr). The Galactic center is in the center of the plot, longitude is increasing from the center to the left. The asymmetry along 
the Galactic longitude derives from the peculiar motion of the sun. 



H — r-Q cos(fe) sm(l)u r — tq sin(fo) sia(£)Ub — (r cos(b) — cos(£)) ut + r cos(fe) cos(^)«0 + (tq — r cos(fe) cos(£)) vq 
In the case of a separable potential given by 

R 2a V£Ro 2a 



ip(R, z) = ip R (R) + tpz(z) , where ip R {R) 



z 2 + D 2 - D 



where G is the universal gravity constant, while Eq, p c s and D are constants, the energies E z and Er obey 

(sin(fe) U r + COs(fe) Uh + Wq) 2 



+ PeSZ 



E z = 



Er = 



(cos(fe) u r — sin(6) M;,) 2 + U 2 + Uq 2 + Vq' 



+ ipa(R) — sin(£) (ui uq — (cos(fe) u r — sin(6) U 6 ) Vq) 



(29) 



(30) 



(31) 



(32) 



(33) 



+ cos(£) ((cos(6) u r — sin(6) Ub) uq + ui vq) , 
while R and z are given by 

R = yj t'q — Itqt cos(fe) cos(£) + r 2 cos(fo) 2 and z = r sin(6) . 

Note that the integration over u r in Eq. (Q) must now be carried numerically since po, fi, K, <t\\ , a z are all functions of u r via Eq. (j29|). 

This model, based on the epicyclic theory, accounts for density and velocity dispersion gradients and is therefore more realistic than 
the Schwarzschild ellipsoid model presented in Subsection [2.2^ The density distribution together with the distribution of the maximum of 
the proper motion along the i coordinate are presented figure projected onto the sphere. The asymmetry along the Galactic longitude is 
produced by the solar motion. 



3 SIMULATIONS 
3.1 Method 

We have chosen to implement a non parametric inversion technique to invert Eq. (^) or (^). The non-parametric inversion problem is 
concerned with finding the best solution to Eq. (Q) or (||) for the underlying luminosity function indexed by kinematic temperature when 
only discrete and noisy measurements of A(b, pb, L) are available (e.g. Dejonghe 1993, Merritt, 1996, Pichon & Thiebaut 1998, Lucy 1994, 
Fadda 1998 and references therein) and most importantly when we have little prejudice regarding what the underlying luminosity function 
should be. In short, the non-parametric inversion corresponds to model fitting in a regime where we do not want to impose (say via stellar 
evolution tracks) what the appropriate parametrization of the model is. It aims at finding the best compromise between noise and bias; in 
effect it correlates the parameters so as to provide the smoothest solution amongst all possible solutions compatible with a given likelyhood. 

An optimal approach should involve a maximum likelihood solution parametrised in terms of the underlying 6-dimensional distribution. 
In practice such an approach turns out to be vastly too costly for data sets involving 10 6 measurements. Binning is therefore applied to our 
ensemble of (£, b, pi,pb, L,B-V) measurements. 
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3.1.1 Non-parametric inversion 

The non parametric solutions of Eq. (j^) and ( |l4| ) are then described by their projection onto a complete basis of p x p functions 

{e k (O e i(P)}k=i,..., P i=i,..., P 

of finite (asymptotically zero) support, which could be cubic B-splines (i.e. the unique C 2 function which is defined to be a cubic over 4 
adjacent intervals and zero outside, with the extra property that it integrates to unity over that interval) or Gaussians: 

v v 
fc=l 1=1 

The parameters to fit are the weights $fej. Calling x = {$ki}k=i,.. P ,i=i,.. P (the parameters) and y = {L 2 cos(b)A(xi, yj)}i=i,.. n> j=i,..n 
(the n x n measurements with L 2 cos(/3) a function of Xi, t/i via Eq. alb), Eq. (ph then becomes formally 



y = ax, 

where a is a (n, n) x (p,p) matrix with entries given by 



Cti,j,k,l 




(35) 



(36) 



i,j,k,l 



For the epicyclic model the measurements are y = {Aijklm = A(£i,bj,fi£ k ^b l ,L m )}i = i,.. ni j = i l ..„ 2 , t k=i,..n 3 J=i,..ni,m=i,..n 5 and 
a is a (ni , ri2 , 7i3 , n4 , 715 ) x (pi , P2 ) matrix with entries given by 



ai,j,k,l,m,q,s = ^ I I I e q (L m r )e 3 (l3)fp(£i,b j , fie k fi bl ,r,u r )r drdu r d/3 > , (37) 

i,j,k,l,m,q,s 

with fp given by Eq. (j26|). 

Assuming that we have access to discrete measurements of Ay ( or Aijkim via binning as discussed above), and that the noise in A can 
be considered to be Normal, we can estimate the error between the measured star counts and the non-parametric model by 



L(x) = x 2 (x) = (y - a-x) x -W-(y - a 



(38) 



where the weight matrix W is the inverse of the covariance matrix of the data (which is diagonal for uncorrelated noise with diagonal 
elements equal to one over the data variance). 

The decomposition in Eq. (^4j) typically involves many more parameters than constraints, such that each parameter controls the shape 
of the function only locally. The inversion problem corresponding to the minimization of Eq. ( p8[ ) is known to be ill-conditioned: Poisson 
noise induced by the very finite sample of stars may produce drastically different solutions since these solutions are dominated by artefacts 
due to the amplification of noise. Some trade-off must therefore be found between the level of smoothness imposed on the solution in order 
to deal with these artefacts on the one hand, and the level of fluctuations consistent with the amount of information in the data set on the other 
hand. Finding such a balance is called the "regularization" of the inversion problem and in effect implies that between 2 solutions yielding 
equivalent likelihood, the smoothest is chosen. In short, the solution of Eq. (D5J) is found by minimizing the quantity 

Q(x) = L(x) + Ai?(x) 

where L(x) and -R(x) are, respectively, the likelihood and regularization terms given by Eq. ( |3ij| ) and 

R(x) = x ± -K-x, (39) 

where K is a positive definite matrix, which is chosen so that R in Eq. ( |39| ) should be non-zero when x is strongly varying as a function of 
its indices. In practice, we use here 



K = K 3 ® I + I ® K 3 + 2K 2 



K 2 



where ® stands for the outer product, I is the identity and K2 = D2 1 ■ D2 , K3 = D3 1 ■ D3. Here D2 and D3 are finite difference second 
order operators (of dimension (p — 2) x p and (p — 3) x p respectively) defined by 



D 2 =Diag 2 [-l,2,-l] = 



-12-10 
0-12-10 
0-12-1 
0-12 



D 3 = Diag 3 [l,-3 ! 3,-l] = 



1 


-3 


3 


-1 








1 


-3 


3 


-1 








1 


-3 


3 











1 


-3 



This choice corresponds a quadratic operator whose kernel include planes and paraboloids. The operator K is typically non zero (and 
therefore penalizes the minimization of Q(x)) for unsmooth solutions (i.e. those leading to strong variations in the coefficients <&ki)- 

The Lagrange multiplier A > allows us to tune the level of regularization. The introduction of the Lagrange multiplier A is formally 
justified by the fact that we want to minimize Q(x), subject to the constraint that L(x) should be in the range iVdata ± v / 27VdItT). In practice, 
the minimum of 



,(40) 



<2(x) = (y-a-x) -W-(y-a-x)+Ax -K-x 



(41) 
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is: 

x = (a ± -W-a + AK) _1 -a x -W-y. (42) 

The last remaining issue involves setting the level of regularization. The so-called cross-validation method (Wahba 1990) adjusts the 
value of A so as to minimize residuals between the data and the prediction derived from the data. Let us define 

a(A) =a-(a ± -W-a + AK) _1 -a ± -W. (43) 

We make use of the value for A given by Generalized Cross Validation (GCV) (Wahba & Wendelberger 1979) estimator corresponding 
to the minimum of 

A () = GCV(\) = min A ( r ll(1 ~ a) " y J|l j • (44) 
[ [trace(l -a)] 2 J 

Note that the model Eq. ( |35| ) is linear and so is Eq. (^) but this need not be the case when positivity is required. We would then resort 
to non linear minimization of Eq. (Mil). 



3.1.2 Positivity 



When dealing with noisy datasets, the non-parametric inversion technique presented above (section 3.1.1) may produce negative coefficients 
in the reconstructed luminosity function. In order to avoid such effects, positivity can be imposed on those coefficients <£>fc; in Eq. (p4|). A 
simple way to achieve positivity is to use an exponential transform and introduce tp so that : 

$ = $ exp (ip) , (45) 
where <l?o corresponds to our first guess for $ (here $o EE 10 3 ). A first order Taylor expansion of Eq. (p^|), together with Eq. (^) yields 
y' ee y - a-<6 = a-$ -<P = ax' , (46) 
which defines y' and x'. We first invert Eq. ( fio] ) for x'. The algorithm is then iterative and we invert in turn for x^ 
y' n = a-x^, where y' n = y - a-<6 n -i and x'„ = $»_i -ip n 
the luminosity function is expressed as 

<$>n = <E>n-l exp (l) CO nvXl/$„-l) (47) 

in Eq. ( fl6| ) for the iteration number n. In practice convergence is controlled via a parameter, ^ conv 6 [0, 1], which fixes the amplitude of the 
correction in Eq. ( ^) in order to remain within the regime of the Taylor expansion. It should be emphasized that using equation Eq. (Q) 
together with Eq. (43) (replacing x by x') does not lead directly to the expected luminosity function but to a correction that has to be applied 
to <E>o . 

We will now proceed to invert Eq. ( Jlo[ ) in two regimes: the Schwarzschild model described by Eq. ( |l8| ) and the epicyclic model, given 
by Eq. ([26]). The former model is dimensionally less demanding, while the latter is more realistic since it accounts for density and velocity 
dipsertion gradients. 



3.2 Simulated Schwarzschild models 

We will first focus on the inversion of Eq. (^), rather than ( A5) or (^fj) (which were shown to be equivalent in the zero asymmetric drift 
approximation) and ( |24| ) (which was also shown to be of the same form). Special emphasis is put on the toy model described in section 2. 1 
while carrying the inversion on a superposition of 4 kinematically decoupled populations with distinct main sequence turn-off magnitudes. 
These are illustrated on Fig. ^ which displays the 4 fictious tracks corresponding to increasing kinematic temperature weighted by some 
IMF on each track. The image in the observed plane (pb, b, L,B-V) of these tracks is shown in Fig. (Q) which shows isocontours of A* 
defined by (corresponding to Eq. with c = 3/4) 

a*,l Afu m ™ /Ssin"(&) [M b 2 sec 2 (&) + 2asin 2 (6)]\ 3/4 

A (b, Hb, L) — A(b, /i b , L) cos(&) 2£"/2-g/3 J ( } 

in the b, fib plane for increasing B-V at a fixed apparent magnitude L — 1/10. The multiple kinematic components of the redder sections 
display distinct extrema for opposite values of yLb at fixed Galactic latitude, b, and also as a function of b at fixed proper motions. In all figures 
7 is chosen equal to 1 (unless specified otherwise) and v to 3. For simplicity, we also numerically approximate K§ in Eq. ( fll| ) by a Gaussian 
since the matrix elements in Eq. (|37j) are then analytic. 




f 

a 



hi S hp 




lowp 



£cJour 

Figure 3. Left: Fictitious tracks corresponding to increasing (from left to right) kinematic temperature Right: decomposition of corresponding colour magni- 
tude diagram into its four components, weighted by the IMF on each track. The image in the observed plane (pt , b, L,B-V) of these tracks is shown in Fig. (W). 



Distribution function 



Potential 



Solar motion 



R, 

R a 
Ro 



p — 2.5 kpc 
2.5 kpc 
= 10 kpc 
= 5 kpc 



o-r q = 48 km/s 
er^ = 24 km/s 
p Q = 0.081 M /pc 3 



D = 240 pc 
E = 48 M Q /pc 2 
p cB = 0.0105 Aio/pc 3 
a = -0.1 



Rq = 8.5 kpc 
V LS R = 220 km/s 
Uq = 9. km/s 
V© = 5.2 km/s 
Wq = 7. km/s 



Table 1. Parameters used for the epicyclic model described in Section E.3 



3.3 Simulated epicyclic models 

In order to test the inversion procedure, a set of 4 HR diagrams with different turn-off luminosity was constructed assuming a mass-luminosity 
relation (MLR) and a Salpeter initial mass function (IMF). The luminosity function of each population scales like 



$o oc L 



(49) 



where r (the slope of the MLR on a logarithmic scale) was set to 3.2, which is caracteristic of the main sequence, and ? to 2.35 (the 
IMF slope). The scaling factor, fixes the number of stars in the simulated galaxy. The tracks associated with those HR diagrams were then 
binned on a 20 x 20 x 4 grid in the [Lq,B — V, beta] space; those HR diagrams represent the absolute luminosity function, $b-v(£o, P). 
The observed counts were then computed assuming that each track corresponds to a given kinematic index and that its distribution can be 
reproduced by the epicyclic model of the same kinematic index i.e.: 



dN(li,bj,m t k,H},,i,L m ,B - V) = o»,i,iM,m, g ,« , ($o(-B - 



!cosbdbdLd(B - V), 



(50) 



where a,i,j,k t i, m ,q,s is given by Eq. Poisson noise was introduced in corresponding histograms used as input for the inversion procedure. 
It should be emphasized that constructing such HR diagrams does not challenge the relevance of our physical model, Eq. (^o|), but only our 
ability to recover a given luminosity function. The model LF need not be very realistic at this stage. The parameters of the epicyclic model 
given in Table [l] were set so as to reproduce the local neighbourhood according to Bienayme & Sechaud (1997) and Vergely et al. (2001). 
Figure |i| show the assumed and reconstruced HR diagrams for the four populations in the [Lo,B — V] plane for this model while figure ^ 
shows the reconstruction error in % for those two figures. 



4 RESULTS 
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Figure 4. A*(b, m^B-V) m the b, /if, plane for increasing B-V (from left to right and top to bottom) at a fixed apparent magnitude L = 1/10 of the model 
described in Fig. (H). Interestingly, the multiple kinematic components of the redder sections display distinct extrema for opposite values of fn, at fixed Galactic 
latitude, b, and also as a function of b at fixed proper motions. 



4.1 the Schwarzschild models 

The above non-parametric inversion technique was implemented on 19 x 19 x 19 data sets (and up to 41 x 41 x 41) corresponding to 
measurements in X, Y, B — V (Eq. (p3|)). For each B — V section, we recover 19 x 19 (resp. 41 x 41) coefficients Xjj corresponding 
to values of U, B, which implies that our resolution in kinematic dispersion is logarithmic. Fig. (Q) shows isocontours of the assumed and 
reconstructed HR diagram as its decomposition in kinematic dispersion. In this zero noise, no bias regime, the relative discrepancy between 
the data and the projection of the model is less than one part in 10 3 while that between the model and the inversion is lower than 10% 
(the corresponding loss in accuracy is characteristic of non-parametric deconvolution). Note that the wiggly structures are a property of the 
model, and are well recovered by the inversion procedure. Fig. (^) shows the actual deprojection overlaid on top of the expected contour of 
the model in the (logarithmic) (/3, L) plane for increasing values of B — V (the projection of the fit in data space is not displayed because 
residuals of the fit would be too small to be seen). Errors in the deprojection are largest for lower contours. Note that the contours of Fig. (Q) 
correspond to sections of the cube shown in Fig. (n) which are orthogonal to those displayed in Fig. (ph. 



4.1.1 Errors in measurements and finite sample 

The above results were achieved assuming infinite numbers of stars and no truncation in apparent magnitude. The Poisson noise induced by 
the finite number of stars (for which accurate photometric and kinematic data is available), as well as the actual error in those measurements, 
are likely to make the inversion of Eq. ^ troublesome. 

Fig. (Q) shows how the error in the recovered HR diagram decreases as a function of the signal to noise ratio in the data which for the 
sake of simplicity, was assumed to be constant while the noise was taken to be Gaussian (corresponding to the large number of stars per bin). 
Note that in reality the signal to noise ratio will clearly be apparent-magnitude dependent, and distance dependent (because of extinction and 
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Figure 5. Left; assumed and Right: reconstructed HR diagram together with its decomposition in kinematic temperature. Note that the wiggly structures 
are a property of the model, and are well recovered by the inversion procedure. The plain dash, dot dash, short dash curves correspond to the 4 dispersions 
associated with the 4 populations with distinct main sequence turn-off radii shown in Fig. (BJ). 



proper motion errors). Fig. (|7|) also shows how the truncation in apparent magnitude induces a truncation in absolute magnitude (here we 
truncate in Y since a truncation in L induces a truncation in Y but none in X given Eq. (p3[)). 

4.2 the epicyclic models 

The inversion technique has been implemented over a36x9x7x7xl0x20x4 model which correspond to a bin size projected onto 
the sphere of 10 x 10 degrees in position sampled linearly, 7 bins in proper motion ranging from -0.2 to 0.2 mas/yr and 20 bins in apparent 
and absolute luminosity correpsonding to an integration over the line of sigth from 0. 1 pc to 4 kpc (those are also linear bins in luminosity, 
which correspond to a logarithmic binning in radius). The four kinematic indexes (ranging from 0.8 to 120) were set to reproduce a series 
of disks with density scale heights ranging from nearly 200 pc to 1 kpc(i.e. corresponding to thin and thick disks). The mean SNR for these 
simulation is 2000, ranging from 20 (on the giant branches) to 70000 on the bottom of the main sequence. 

Figure ^ shows the reconstruction error in the [B — V, Lq] plane correpsonding to the HR diagram shown in Fig [| The main sequence 
and the different turn-off are well reconstructed (the error lies well below 1% for the faint part of the main sequence and is less than 10% at 
the turn-offs). The Red Giant Branches (RGB) are also well reproduced even though it strongly depends on the age of the population (via /?). 
This can be understood if we look at the number of stars in the different regions on the [B — V, Lq] plane. Older ( resp. younger) populations 
have larger (resp.lower) number of stars on the RGB and the SNR is increasing (resp. decreasing) correspondingly. We note that the four 
tracks are recovered without creating any spurious structure. The luminosity functions (Lq) are recovered within 1% uncertainty (in mean 
value) for the oldest population and within 20% for the youngest (note that sometimes, the reconstruction error increases up to 100% when 
no stars are recovered on the RGB) 



5 DISCUSSION AND CONCLUSION 

The main result of this paper is a demonstration that the generalized stellar statistic equation including proper motions, Eq. (^|), can be 
inverted, giving access to both the kinematics and the luminosity function. The inversion was carried for two rather specific functional 
decompositions of the underlying distribution (namely, constant ratio and possibly singular Schwarzschild ellipsoids plane parallel models) 
and a more realistic physical model (the epicyclic Shu model) which accounts for gradients. The inversion assumes that the departure from 
harmonicity of the vertical potential, and/or the asymmetric drift or the Sun's vertical velocity, wq, are known. Indeed the break in the 
potential yields a scale which reflects the fact that statistically the dynamics (i.e. the velocities) gives a precise indication of distances in 
units of that scale. The asymmetric drift or vertical component of the sun's velocity provide another energy scale (and therefore a distance 
scale). The existence of more than one distance scale is mathematically redundant, but practically of interest for the purpose of accounting 
for local and remote stars. 

In a nutshell, it was shown in section ^ that Eq. (Q) has solutions for families of distribution obeying Eq. (Q) (singular ellipsoid) or 
Eq. (|l8|) (Schwarzschild ellipsoid). Those solutions are unique and can be made explicit for a number of particular cases: Eq. (JlTp (pin-like 
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Figure 6. assumed model (plain line, filled contour) and non parametric deprojection (dashed line) overlaid on top of the expected contours in the 
(log 1/L, log/3) plane for increasing B — V sections. Errors in the deprojection are largest for lower contours. Note that these sections are orthogonal 
to those superposed in Fig. (|^). 



velocity ellipsoid), Eq. (A5) (constant ratio /3r/ (3 z , Wq ~ 0), Eq. (48) (constant ratio /3r//3 z and P^/Pz either with vq ~ v$ ,wq ~ 0, 
or vq — 7^ 0,uiq 7^ & x ~ 0: statistical sec ular p arallaxes). In all other instances, the solution can be found via the general 

the only constraint being the computation of the model matrix generalizing 



3.1.1 



non- par ametric inversion procedure described in section 
Eq. (p7p (which might require numerical integration, as shown for instance in section 2.3); in this more general framework it remains also 
to demonstrate that the inversion will converge towards a solution which is unique. For instance, in the regime where the epicyclic model 
has been tested (subsection |3.3| ) a unique solution seems to be well defined. The luminosity function of each kinematical component is well 
recovered throughout the HR diagram. 

More tests are required before applying the method to real data, and are postponed to a companion paper (Siebert et al. in preparation). 
For a given vertical potential, it appears that the modelling of star counts indexed by proper motion A\(m, [ii, [ib',l,b) has a solution 
for most model parameters. Many different models based on distinct priors have produced realistic magnitude and colour star counts but 
failed to predict proper motion measurements accurately (for instance, note that the Besancon model - which relies on a nearly dynamical 
consistent model - produces a good fit to proper motion surveys (Ojha, 1994), while dynamically inconsistent models are more problematic 
(Ratnatunga, 1989)). 



It should be emphasized that the inversion method presented in section [3.1.l| is a true deconvolution and should give access to a kinemat- 
ically indexed HR diagram. Together with some model of the time evolution of the different kinematic components (via, say, a disk heating 
mechanism), the indexing could be translated into one on cosmological time, hence providing a non-parametric measurement of the local 
neighbourhood luminosity function which is complementary to that obtained by evolutionary track fitting with an assumed initial mass func- 
tion and star formation rate (see, for instance, Hernandez, Valls-Gabaud & Gilmore, 1999). Note that, conversely, the agreement between the 
standard direct method to predict the local luminosity function and the method presented here could be used to measure the Galactic potential. 



The deepest photometric and proper motion of whole sky survey available is the Tycho-2 catalogue ( H0g et al., 2000) which is a new 
reduction of the Tycho data (H0g et al., 1998). Many Tycho stars are disk giants and subgiants covering a large range of distances; the 
method developed here can be applied to these stars and will allow us to recover their luminosity function without any prior information 
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Figure 7. Left Panel: the mean absolute residual of the luminosity function, |$™oov _ <E>"^ put |/ h^;" 10 "* versus the signal to noise ratio in 

logarithmic coordinates. This graph demonstrates that the non-parametric inversion sketched in section 3.1.1 is robust with respect to sampling or measurement 



noise. Rig/it Panel: the effect of truncation in magnitude on the main sequence: plain line: recovered HR diagram with a truncated data set; dashed line: 
recovered HR diagram without truncation. As expected, the truncation in apparent magnitude removes the information at the bottom of the main sequence. 
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Figure 8. Left: Assumed HR diagram for the epicyclic model. The four populations have distinct turnoff point and kinematic index. Right: Reconstructed HR 
diagram. Lq is expressed in unit of Lq. Note that all four populations are well recovered. Main sequence and turn-off are reconstructed within 10% error (less 
than 1% for the lower part of the main sequence due to the large number of stars in that part of the HR diagram). Giant branch are also recovered though the 
reconstruction error is higher. 



from stellar evolution tracks. We intend in a forthcoming paper to apply the method presented here to the Tycho-2 catalogue (H0g et at., 
1998) and to other proper motion catalogues in order to determine the luminosity function of stars in the Solar neighbourhood. We will 
investigate the limitations introduced by a magnitude-limited catalogue, by the finite size of catalogues and also by our limited knowledge 
of the Galactic potential. Reddening is also bound to be a concern since it will bias apparent luminosities as a function of i and b. If the 
reddening is diffuse and the absorbing component law is known, the kernel of e.g. Eq. (|l0 1 will simply be modified accordingly. Alternatively, 
multicolour photometry could be sufficient to constrain the spatial extinction law. Of course the dimensionality of the problem is increased 
by the number of different colour bands used, since the analysis must be carried while accounting for all colours simultaneously. 

The final error on the recovered luminosity functions will depend on the photometric errors of the observational catalogued 0.1 for 
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Figure 9. Left panel: Reconstrution error in the [B — V, Lq] plane for the epicyclic model shown in Figs. ^ Right panel: Model vs. recovered luminosity 
function &p(Lo) for the four kinematic indexes correponding to the oldest population (lower curve) to the youngest (upper curve). The LFs are plotted on 
a logarithmic scale and arbitrary normalised. The curves corresponding to the two kinematic index where shifted along the y axis. Plain lines correspond to 
the model LF while dash-dotted lines are the reconstructed LF. Note that the LF is well reconstructed for the main sequence (at low luminosity) and for the 
turn-off. The total luminosity function summed over the kinematic index is also displayed as the top thick lines. The bumps at low and high luminosity are 
properties of the model and correspond to the lower part of the main sequence and to the subgiant branch of each population. 



the Tycho-2 catalogue down to 0.013 for Vt < 9, ~ 0.05 — 0.10 for photographic surveys). It will also depend on the relative proper 
motion accuracy (S(m) ~ 25(/i)/u M , with cr M the typical dispersion for a stellar group at a given distances). With the Tycho-2 catalogue 
completed by proper motions (with an accuracy of 2.5 mas/y), and for disk giants with velocity dispersions from 10 to 50 km/s and proper 
motion dispersions from 2 to 10 mas/y, the accuracy on the recovered luminosity function will be limited to about 0.5 magnitude. Closer 
(and fainter) stars with proper motions from photographic catalogues will constrain the lower part of the luminosity function with a higher 
accuracy. 

In a decade, sky surveys by the Fame, Diva and GAIA satellites will probe the Galactic structure in superb detail, giving directly access 
to larger volumes of the 6D stellar phase space of the Galaxy. It will remain that farther out, only proper motions and photometry will have 
sufficient accuracy and generalization of methods such as that derived here will be used to extrapolate our knowledge of the kinematic and 
luminosity function of the Galaxy and its satellites. For instance Appendix [b| sketches the possible inversion of an external globular cluster 
luminosity function with GAIA quality photometry. 
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APPENDIX A: ASYMPTOTIC ANALYTIC SOLUTIONS FOR THE SCHWARZSCHILD MODEL 



Let us demonstrate that Eq. (0) has explicit analytic solutions for families of distribution obeying qlq), using the inversion procedure sketched 



in section 2. 1 



Al Slices towards the Galactic centre 

For the sake of simplicity let us first restrict the analysis to uq — vq — wq — and assume first we have measurements only in the direction 
I = 0. The integration over u r then yields 



f (r,u)du r = y ^pj- exp 

tq — r cos(6) 



where 

/3- 1 =/3 fl 1 sin 2 (6) + 



R 



/3 Z 1 cos 2 (6) . 



Without loss of generality let us integrate over ui 

Pi, 



ffj(r, u)du r due 



■ exp 



{~-\PbU 2 b -f3 z il) z (z) 



At large distances to the Galactic centre, both R and r© are large compared to r and Eq. ( A2 ) becomes 

ft" 1 =^ 1 sin 2 (6)+/3- 1 cos 2 (6). 



(Al) 



(A2) 



(A3) 



(A4) 



Let us now also assume that Pr and f3 z are known monotonic functions of a unique parameter (3. We may now convolve Eq. (A3) with 
the luminosity function sought, $ \Lr 2 , /?] so that 



A[b, fi b , L] 



[ir 2 ,/?] exp ^-ir- 2 /3(,/i 2 ~ (3 z ip z (rsm(b))j r 3 drd/3 . 



(A5) 



Eq. (Q) appears now as a special case of Eqs. (A4)-(A5l corresponding to (3r — > 00. Even though the convolution in Eq. ( ^ ) is less 
straightforward than that of Eq. (^|), and so long as ip z is not purely harmonic, Eq. ( A5) will have a non-trivial solution for In particular, 
if the ratio of velocity dispersions (3r/ f3 z is assumed constant, Eq. (H) still holds but with f3 — [3 Z , and x replaced by x' defined by 



sin 2 (6) 
x - a -i-i + 



2Lcos 2 (fo) +C2Lsin 2 (6) 



where ^ 



with A'[b,ii b ,L] = A[b, fi b , L]^l + ^ta,n 2 (b) . 



Note that if vq and wq are not negligible, Eq. (|A5|) becomes 



A [b, fi b , L] 



[Lr 2 ,/3] exp (--^/3b(ub + cos(b)w e ) 2 - /3 z tp z (r sin(fo))^ r 3 drd(3 , 



(A6) 



and is of the form discussed below as Eq. ( A8 ) with i = 



A2 Slices away from the Galactic centre 



For any direction t 7^ when tq — ■» 00, the kinetic dispersion (replacing in Eq. ( A4 )) along Galactic latitude is then given by 

(i b x = (Pr 1 cos 2 t + f)^, 1 sin 2 1} sin 2 b + P^ 1 cos 2 b , 
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and Eq. (ch is replaced by 



A[b,£,fib,L] = j J y ^-c£> [Lr 2 ,/?] exp {—(3b{rpb +cos(fo)™ — sin(6) sin(Z) [v@ — v^]) 2 f3 z ~ip z (r sin(6))) r 3 drd/3, 

which can be rearranged as (again with (3 = (3 Z ) 

7 



L cos{b)A 2 \b, £, Hb, L] 

with z 2 given by Eq. ( [25^ 

sin 2 (6) 



$ [u 2 , 0\ exp {-Pu 2 x 2 + /3uy 2 ~ /3z 2 - Px{uy)) u 3 dudf3 , 



x 2 = a- 



2/2 = 



£ 2L cos 2 (6) + 2L sin 2 (&) (£ R cos 2 (£) + ^ sin 2 (I)) 
Hb ( [v© — v<j>] sin fe sin I — wq cos 6) 



v/Zcos 2 (&) + VI sin 2 (6) (£flCos 2 (£) + ^sin 2 (£)) ' 
A a [M,W.,£] = A[b, fj, b , L] y/i + tan 2 (&)[^ cos 2 ^) + ^ sin 2 ^)] . 



where £ fl = — , ^ = — . 

PR P<t> 



(A7) 

(A8) 

(A9) 
(A10) 
(All) 



In the region where the asymmetric drift and the z-component of the Sun's velocity can be neglected, vq w v<p and Wq « 0, j/2 
and ^2 vanish and Eq. ( A8) is formally identical to Eq. (JToj) ; once again the solution of Eq. ( A8 ) is g iven by Eq. ( |T^ ) with the appropriate 
substitutions. Alternatively, in the regions where either wq or vq — cannot be neglected, Eq. (A8) has a unique solution even if \ = 0, 
which can be found along the section pb = (Note that when — > oo, we can always a ssume uq = by changing the origin of Galactic 
longitude, £). Indeed, Eq. (48) becomes Eq. ( p^j ) which is of the form described in section 2.1.2 with v — 0, X3 replacing x and z 2 replacing 
y; the corresponding solution is found following the same route. It is analogous to statistical secular parallaxes (note nonetheless that the 
section fib = might not be sufficient to carry the inversion without any truncation bias since log (22) spans ] — 00, Z[ when b and £ vary 
with Z a function of £r,£*, wq and vq — v$). 

Turning back to Eq. (48), it remains that for more general \ it can still t> e inverted via the kernel, K 2 (x 2 , y 2 , z 2 ,y\u, (3), which depends 
explicitly on \' 



K 2 (x 2 ,y 2 ,z 2 ,y\u,(3) 



exp {-(3ux 2 + (3uy 2 - f3z 2 - f3x{uy)) u 3 



Note that the multi-dimensionality of the kernel, K 2 , is not a problem from the point of view of a x non-parametric minimisation described 
in section H. 



APPENDIX B: EXTERNAL SPHERICAL ISOTROPIC CLUSTERS 

Consider a satellite of our Galaxy assumed to be well-described as a spherical isotropic cluster with a luminosity function indexed by this 
kinematic temperature. Let 4-7T 2 A\(pb, L, R)pdpdLRdR be the number of stars which have proper motions, p 2 — pi + p 2 , apparent 
luminosity L at radius R from the centre at the wavelength A. This quantity is a convolution of the distribution function f(e, (3) (a function 
of energy, e, and (3 = 1/cr 2 ) and the luminosity function, g\(Lo,(3), a function of the intrinsic luminosity, Lo, the population, (3, and 
wavelength A: 

A x (im,L,R) = J J J f(e,l3)g x (f3,Lr' 2 )df3dzdv z , (Bl) 
which can be rearranged as 

A x (f,,L,R)=4 [([ f(s,P)g x (f3,Lr' 2 ) J^_ —= d£ df3 , (B2) 

J J J Vr 2 -R 2 y /2(^ + e)-v 2 

where v 2 — p 2 r' 2 is the velocity in the plane of the sky, and r' 2 — (r 2 — R 2 + t-q) the distance to the observer and tq the distance to 
the cluster, and r the distance to the cluster centre. The potential can be derived non-parametrically from the projected density (using Jeans' 
equation). Indeed the enclosed mass within a sphere of radius r reads 

M dyn «r)=r 2 ^ = - r -^fl, (B3) 
dr p dr 

where ip(r) is the gravitational potential p(r) the density and <r(r) the radial velocity dispersion. The surface density is related to the density 
via an Abel transform: 

00 

/* /*°° rl 

E(R) = j p(r)dz = 2 j p(r) = A R (p) , (B4) 
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where E(-R) is the projected surface density and R the projected radius as measured on the sky. Similarly the projected velocity dispersion 
<Tp is related to the intrinsic velocity dispersion, a 2 (r), via the same Abel transform (or projection) 



E(R)<4(R) = 2 / P (rV(r) 



rdr 



Vr 2 - R 2 



A R (pcr 2 



(B5) 



Note that E(.R)<7p is the projected kinetic energy density divided by three (corresponding to one degree of freedom) and p(r)a 2 the kinetic 



(B6) 



energy density divided by three. Inserting Eqs. (B4 )-(|B5|) into Eq. (B3 1 yields 



M dyn (< r) = 



&Ar \Val) ^ wh . le p{r) = l d M ^ (< r) ^ = _^ Qp 



Ar 1 ^) dr ' ry ' ' Anr 2 dr 

The underlying isotropic distribution is given by an inverse Abel from the density. 



fie) 



d 2 p dip 



VEtv 2 J dtp 2 Ve^ip 



F(p)exp(-pe)dp 



(B7) 



where an isothermal decomposition over temperature j3 was assumed for the distribution function (this assumption is not required: any 
parametrized decomposition is acceptable). So 



FQ3) = C' 1 [/(e)] = C- 1 [A' 1 (A' 1 (E))] 

where C is the Laplace operator. 
Calling 



(B8) 



G[Y] = j^-*^ =V^EMVY)e- Y , g^,Lr' 2 )=g,{f3,Lr' 2 )F{P)p- 3/2 



(B9) 



Eq. (B2h becomes 



A x {fJ,,L,R) = 2V2 J (J G 



ff i [p,L(r 2 -R 2 +r%)] 



rdr 



Vr 2 - R 2 



dp 



(BIO) 



where G is a known kernel while g\ is the unknown sought luminosity function. Eq. (BIO) is the direct analogue to Eq. (Q). It will be 
invertible following the same route with GAIA photometry. (With today's accuracy in photometry, for a typical globular cluster at a distance, 
Tq of, say 10 kpc, the relative positions within the cluster are negligible w.r.t r Q : r 2 — R 2 <C , therefore 



A x (v,L,R) = 2V2j (J G 



rdr 



Vr 2 - R 2 



9l [p,Lrl] dp 



L is then also mute, and the inversion problem shrinks to one involving finding the relative weights, ql [P] of a known distribution). 



